ON SCALE FUNCTIONS OF SPECTRALLY NEGATIVE LEVY PROCESSES WITH 

PHASE-TYPE JUMPS 

MASAHIKO EGAMI AND KAZUTOSHI YAMAZAKI 

Abstract. We study the scale function of the spectrally negative phase-type Levy process. Its scale 

function admits an analytical expression and so do a number of its fluctuation identities. Motivated by the 

fact that the class of phase-type distributions is dense in the class of all positive-valued distributions, we 

^~~^ propose a new approach to approximating the scale function and the associated fluctuation identities for a 

^^ general spectrally negative Levy process. Numerical examples are provided to illustrate the effectiveness 

CnJ of the approximation method. The extension to the case with jumps of infinite activity is also discussed. 
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Q In the last decade, significant progress has been made regarding spectrally negative Levy processes, and this 

is mainly due to the scale function. As can be seen in the work of, for example, [7, 25], a number of fluctuation 

'^ identities concerning spectrally negative Levy processes can be expressed in terms of scale functions. There are 

■^ numerous applications in applied probability including optimal stopping, queuing, branching processes, insurance 

^ and credit risk. Despite these advances, a major obstacle still remains in putting these in practice because scale 

O functions are in general known only up to their Laplace transforms, and only a few cases admit explicit expressions. 

lO Typically, one needs to rely on numerical Laplace inversion in order to approximate the scale function; see [24, 38]. 

Q In this paper, we propose a phase-type (PH)-fitting approach by using the scale function for the class of spectrally 

^—1 negative PH Levy processes, or Levy processes with negative PH-distributed jumps. Consider a continuous-time 

^ Markov chain with some initial distribution and state space consisting of a single absorbing state and a finite number 

S^ of transient states. The PH distribution is the distribution of the time to absorption. The class of PH distributions 

3 includes, for example, the exponential, hyperexponential, Erlang, hyper-Erlang and Coxian distributions; see, 

e.g. Section 3 of [2]. It is known that the class of PH distributions is dense in the class of all positive- valued 

distributions, and consequently the scale function of any spectrally negative Levy process can be approximated by 

those of PH Levy processes. 



1. Introduction 
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Our objective is to obtain the scale function of the spectrally negative PH Levy process, and use it to approximate 
for a general case. We achieve its closed- form expression by taking advantage of the structure of its Wiener-Hopf 
factor as obtained in [3]. In order to verify that this is indeed a useful alternative to the Laplace inversion approach, 
we show convergence results via the continuity theorem and then conduct a series of numerical experiments. 

One major advantage of this approach is that the approximated scale function is given as a function in a closed 
form, which enables one to analytically differentiate/integrate to obtain other fluctuation identities explicitly. The 
derivative of the scale function is used to identify, for example, the fluctuation of its reflected process and is 
commonly applied in insurance risk literature [5, 27, 30]. The integration of the scale function with respect to the 
Levy measure is required when overshoots/undershoots are involved. On the other hand, these calculations must 
be conducted numerically when the scale function is approximated via numerical Laplace inversion; because the 
scale function increases exponentially and its derivative tends to explode near zero, minimizing numerical errors 
tends to be a difficult task. 

Another advantage of the PH- fitting approach is that the Laplace transform of the PH distribution has an expUcit 
expression. The Laplace inversion approach inverts the equality written in terms of the Laplace exponents which 
admit analytical expressions only in special cases (see (2.1)-(2.2) below). In other words, it contains two types 
of errors: 1) the approximation error caused while computing the Laplace exponent and 2) the error caused while 
inverting the Laplace transform. On the other hand, the PH-fitting approach only contains the PH-fitting error 
thanks to the closed-form Laplace transform of the PH distribution. 

The PH-fitting approach is particularly efficient when the Levy density is completely monotone because the 
fitting can be conducted by hyperexponential distributions (i.e. hyperexponential-fitting) and there exist algorithms 
that are guaranteed to converge (see, e.g., [15]). The class of Levy processes with completely monotone Levy den- 
sity is rich. It includes compound-Poisson processes with long-tailed distributed jumps such as the Pareto, WeibuU 
and Gamma distributions and other processes such as variance gamma [35, 36], CGMY [11], generalized hyper- 
bolic [14] and normal inverse Gaussian [6] processes. The completely monotone assumption of the Levy density 
is sometimes required, for example, in the optimal dividend problem [33, 34]. In order to give examples of how 
hyperexponential-fitting can be applied to compute fluctuation identities, we obtain the closed-form expressions of 
the overshoot/undershoot distributions at the first down-crossing time for the hyperexponential case, and use it to 
approximate those for the processes with Weibull/Pareto-type jumps. 

The Levy measure of the PH Levy process is only limited to a finite measure. To complement this, we also 
consider the meromorphic Levy process [22], whose Levy measure may be infinite and is given as a countable 
sum of exponential functions. This generalizes the hyperexponential Levy process as well as a number of Levy 
processes such as Lamperti-stable processes [10, 13], hypergeometric processes [23, 26] and those in the /3- and 0- 
families [20, 21]. The corresponding scale function can be expressed explicitly as a countable sum of exponential 
functions, which can be effectively approximated by a finite sum with some analytical error bounds. We show 
numerically the effectiveness of the approximation procedure using examples with /3-processes [20]. 

The rest of the paper is organized as follows. Section 2 studies the spectrally negative PH Levy process and 
its scale function. We show via the continuity theorem that the PH-fitting approach can approximate the scale 
function for a general spectrally negative Levy process. We further obtain, for the hyperexponential case, the joint 
distribution of the overshoot and undershoot at the first down-crossing time. Section 3 obtains the scale function 
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for the meromorphic Levy process and its upper and lower bounds. Section 4 verifies the effectiveness of the 
PH-fitting approach through a series of numerical results. All proofs are given in the appendix. 

2. SPECTRALLY NEGATIVE PHASE-TYPE LEVY PROCESSES 

2.1. Scale functions. Let (il,J^, P) be a probability space hosting a spectrally negative Levy process X = 
{Xt] t > 0}, F^ be the conditional probability under which Xq = x, and F := {Tt : t > 0} be the filtration 
generated by X. The process X is uniquely characterized by its Laplace exponent 

1 r^ 

(2.1) i){s) ■.= \og¥P [e'^^] = jls + -a'^s^ + I (e"^ - 1 - szl|,>_i})n(dz), sGC 

^ J -oo 

where 11 is a Levy measure with the support (— oo, 0) and satisfies the integrability condition L_^ qn (1 Az^)n(d2;) < 
oo. It has paths of bounded variation if and only if 



0- = and / (1 A |z|)n(d2;) < oo; 



'(-00,0) 

see, for example, [25], Lemma 2.12. In this case, we can rewrite the Laplace exponent (2.1) by 

r-O 



with 



V'(s) = ^s+/ (e'^-l)n(dz), sGC, 

J —oo 

rO 

fi := pL — / zll{dz). 

Fix q > and any spectrally negative Levy process with its Laplace exponent ^. The scale function VF^'^^ 
M I— )• R is a function whose Laplace transform is given by 

/>oo 1 

(2.2) / e-^"H^('?)(x)dx = — , s>Cg 

Jo ^(s) - Q 

where 

(2.3) Cg := sup{s > : i;{s) =q}, q> 0. 

We assume W^'^\x) = on (— oo, 0). 

Let us define the first down- and up-crossing times, respectively, by 

(2.4) T~ := inf {t > : Xt < a} and t^ := inf {t > : Xt > b} , a,b e M. 
Then we have for any x < b 

(2.5) 

e'^'^^l/ +^ - +^ 1 ^ ^ ' ("^j and E^le-'^^oi ] = Z^^\x) - Z'^'^'> (b) ^ ^ ,^'^} 



where 



(2.6) Z(«)(x) :=l + q [ l^('')(y)dy, x G M. 

Jo 

Here, we disregard the case when X is a negative subordinator (or decreasing a.s.). 
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We also consider a version of the scale function W(; : M i— )• M that satisfies 

(2.7) M/(«)(a;) = e^^^VF^,(x), xGM 
with its Laplace transform 

(2.8) r e-''^W^^{x)dx= -^ , s>0. 

Jo i^{s + Cq)-q 

Suppose Pc, for any given c > 0, is the probability measure defined by the Esscher transform 

(2.9) ^ 

^ ' dP 



= e^^t-^i'^)\ t > 0; 



see page 78 of [25]. Then W(^^ under P^^ is equivalent to W^^^ under P. Furthermore, it is known that Wq is 
monotonically increasing and 

(2.10) W^^{x)/^{i,'{Qr^ asx^oo, 
which also implies that the scale function M/w) increases exponentially in x; 

(2.11) W-'^Hx)'^^y^ asx^oo, 

ViCq) 

except for the case q = and ip'{Co) = 0. 

Regarding the smoothness of the scale function, if its jump distribution has no atoms, then Ty('') G C^(0, oo); 
if it has a Gaussian component {a > 0), then W^'^'' G C^i^, cc); see [12]. In particular, a stronger result holds 
for the completely monotone jump case. Recall that a density function / is called completely monotone if all the 
derivatives exist and, for every n > 1, 

(-l)"/(")(x)>0, x>0, 

where f^^> denotes the n*^ derivative of /. 

Lemma 2.1 (Loeffen [33]). If the Levy measure has a completely monotone density, W',-- is also completely mono- 
tone. 

Finally, the behavior in the neighborhood of zero is given as follows; see Lemmas 4.3-4.4 of [29]. 

Lemma 2.2. For every q > 0, we have 



f,A , , 0, unbounded variation / ,w , 

M^('?)(0) = <^ ,' , , ) and W^'^^ (0+) 

I - , bounded variation I 



^, a>0 

oo, (7 = and n(— oo, 0) = oo 

- — 2°°' ' , compound Poisson 
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2.2. Spectrally negative PH Levy processes. Consider a continuous-time Markov chain Y = {Yt; t > 0} with 
finite state space {1, . . . , m} U {A} where 1, . . . , ?n, are transient and A is absorbing. Its initial distribution is 
given by a simplex a = [ai, . . . , am] such that Oj = P {Yq = i} for every i = 1, . . . ,m. The intensity matrix Q 
is partitioned into the m transient states and the absorbing state A, and is given by 

T t 





Q- 



Here T is an m x m-matrix called the PH-generator, and t = —Tl where 1 = [1, . . . , 1]'. A distribution is called 
phase-type (PH) with representation {m, Q, T) if it is the distribution of the absorption time to A in the Markov 
chain described above. It is known that T is non-singular and thus invertible; see [1]. Its distribution and density 
functions are given, respectively, by 

F{z) = 1 - ae^""! and f{z) = ae^H, 2 > 0. 

Let X = {Xt;t > 0} be a spectrally negative Levy process of the form 

Nt 

(2.12) Xt- Xq = ^it + aBt-^Zn, < t < 00, 

n=l 

for some /i G M and cr > 0. Here B = {Bt] t > 0} is a standard Brownian motion, A^ = {Nf, t > 0} is a Poisson 
process with arrival rate A, and Z = {Z„; n = 1, 2, . . .} is an i.i.d. sequence of PH-distributed random variables 
with representation {m, a, T). These processes are assumed mutually independent. Its Laplace exponent is then 

(2.13) Tpis) = HS + -a'^s^ + A [q.{sI - TyH - l) , 

which is analytic for every s e C except for the eigenvalues of T. 

Disregarding the negative subordinator case, we consider the following two cases: 
Case 1: when cr > (i.e. X is of unbounded variation). 
Case 2: when cr = and /i > (i.e. X is a compound Poisson process). 

Here, in Case 2, down-crossing of a threshold can occur only by jumps; see, for example. Chapter III of [7]. On the 
other hand, in Case 1, down-crossing can occur also by creeping downward (by the diffusion components). Due to 
this difference, the form of the scale function differs as we shall see. 
Fix (7 > 0. Consider the Cramer-Lundberg equation 

(2.14) ^P{s)=q, 

and define the set of (the absolute values of) negative roots and the set of poles: 

Zq := {i : i^{-Ci,q) = q and 7^(C^,g) > 0} , 



Jq := \j : ^ = and n{ri,) > o] 



The elements in Xg and Jq may not be distinct, and, in this case, we take each as many times as its multiplicity. By 
Lemma 1 of [3], we have 



IX, 



qi 



\Jq\ + l, for Case 1, 
|j7q|, for Case 2. 
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In particular, if the representation is minimal (see [3]), we have \Jq\ = m. 

Let Cq be an independent exponential random variable with parameter q and denote the running maximum and 
minimum, respectively, by 

Xt = sup Xg and X_^ = inf Xg, t > 0. 

0<s<t 0<s<t 

The Wiener-Hopf factorization states that q/{q — ip{s)) = ip'^{s)(p~{s) for every s G C such that TZ{s) = 0, with 
the Wiener-Hopf factors 

(2.15) v.,(s):=E[exp(sX,j] and (/p+(s) := E [exp(sXeJ] 

that are analytic for s with Tl{s) > and Tl{s) < 0, respectively. By Lemma 1 of [3], we have, for every s such 
that n{s) > 0, 

from which we can obtain the distribution of X^ by the Laplace inverse via partial fraction expansion. 

As in Remark 4 of [3], let n denote the number of different roots in Zg and m,j denote the multiplicity of a root 
^i^q for i = 1, . . . , n. Then we have 

i=l k=l 

where 



n mi ,^ \k—l 

iiAi) F {-x,^ e dx} = ^ J: 4'^^^'^%^^"'"'^'d^' ^ > 



^(fc) ._ 1 d^'-'' V9-(g)(g + e,,,)"' 



*'"? ■ {rrii - k)l ds"'-'' i 



t,q 



Notice that this can be simplified significantly when all the roots in Zq are distinct. 

2.3. Scale functions for spectrally negative PH Levy processes. We obtain the scale function focusing on 
the case q > 0. The scale function when q = and X drifts to +oo can be derived by using W^'^^x) = 
IP"^ {±Koo ^ 0} /V''(0) ^nd the ruin probability (19) of [3] by taking g — )■ 0; Kyprianou and Palmowski [27] briefly 
stated the scale function when q = 0, X drifts to +cxd and all the roots in Zq are distinct. The case q = and 
X drifts to — oo can be obtained indirectly by change of measure (2.9); in this case, W^^\x) = e^''^W^„(x) with 
Co > and X drifts to +oo under P^^. 

Before obtaining the scale function, we shall first represent the positive root Cq (2.3) in terms of the negative 
roots {^i^q] i G Zq}. Let us define 



(2.18) gg:=Y^A^^gC,,q, q>0, 



and by Lemma 2.2 



(2.19, ^--CWiM^Wi'y(0^)4h,^/ZcZ'2 
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Lemma 2.3. For every q> Q, we have 

Q Qq 

We now obtain the version of the scale function W(^ (•). In the lemma below, Wq (0) = W^'^\Q) is either or 
- depending on if it is Case 1 or 2; see Lemma 2.2. 

Lemma 2.4. For every q> Q, we have 



w^c,(-)-^c.(o) = ^EE4? 



^ n rrii , ^ 



1=1 k=\ 



,i,q 



Qq + S,i,q 



k-1 



l_e'(C1+i^,q)^y^ 



iiCq + ^^,q)xy 



3=0 



j! 



x> 0. 



Lemmas 2.2 and 2.4 and (2.7) show the following. By Lemma 2.3, these expressions can also be rewritten in a 
different way. 

Proposition 2.1. Fix q > and x > 0. 
(1) For Case 1, 






^ i=\k=l ^'"^ VCq+Cr,<? 



o<iq^ _ o-i 



k~\ 

j=o 



(.{(q + ^^,q)^y 



(2) For Case 2, 



^ n mi , ^ 



h^W(x) = ^eea^: 



i=lk=l " VCg + C*,g 



k-1 



oCq^ _ ^-6,,a^ "^ \\^q^^t,q 



j=0 



{{(:q + i^,q>y 



j! 



+ -e 



Cqa; 



The scale functions obtained above are infinitely differentiable. In particular, the first derivative becomes 

k-1 '=-2 

';,i,q)->-l 

(k 






i,q 



j=l fc=l 



Cg + Ci,, 



Coe^'^ 



^ ' j=0 



{{Cq+Ci,q)xy 



J! 



for Case 1 and 



^(g)'(^^ - ''^ 






EEa 



i=ifc=i V^g + ?«,g 



Cge'^""' + 6,ge ^'•'' 



(fe-1)! 



fe-2 



Qe-^'-^E 



((Cg + C.,.)^^ 



i=o 



j! 



+ -Cqe^"' 



for Case 2. 

When all the roots in Iq are distinct, the scale functions above can be simplified and have nice properties as 
discussed in the following corollary. 



Corollary 2.1. Suppose all the roots in Xq are distinct. 
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(1) The scale function can be simplified to 

(2.20) 






iein 



for Case 1 and 2, respectively where 



gCg^ g ?i,ij2 



l(l) 



+ -e 



Cq^ 






-'«,9 



9 Cg + ii,q ■ 



(2) 1^(5)' /^ convex. 

(3) W^ /5 completely monotone. 

In the lemma above, by (2.10)-(2.1 1), we must have 

(2.21) 



v^y^S""^^^ ^"' ^ = ^^-+^' 



ieXg 



«ei„ 



^ 



respectively for Case 1 and 2. 



Example 2.1 (Hyperexponential Case). As an important example where all the roots in Xq are distinct, we consider 
the case where Z has a hyperexponential distribution with density function 

m 

/(z) = ^p,7?,e-^^^ z>0, 
for some < ryi < • • • < 7]^ < oo. Its Laplace exponent (2.1) is then 



(2.22) 



^ m 

V'(s) =/is+ ^cr'^s^ -X'^pj — 



^ ^^• + ^ 



Notice in this case that —iji, ... , —ijm are the poles of the Laplace exponent. Furthermore, all the roots in Xq are 
distinct and satisfy the following interlacing condition for every q > 0: 

(1) when o" > 0, there are m+ 1 roots — ^i,g, • • • , —^m+i,q such that 

< 6,9 < ^1 < 6,g < • • • < ??m < Cm+l,q < OO; 

(2) when a = and /U > 0, there are m roots — 6,gi • • • i ~Cm,q such that 

< Cl,g <m < 6,g < • • • < S,m,q < r/^ < OO. 

77je cZfl'i'i' of hyperexponential distributions is important as it is dense in the class of all positive-valued distributions 
with completely monotone densities. 
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2.4. Approximation of the scale function of a general spectrally negative Levy process. The scale function 
obtained in Proposition 2.1 can be used to approximate the scale function of a general spectrally negative Levy 
process. 

By Proposition 1 of [3], there exists, for any spectrally negative Levy process X, a sequence of spectrally 
negative PH Levy processes X^") converging to X in L)[0, oo). This is equivalent to saying that X|" — ;• Xi in 
distribution by [18], Corollary VII 3.6; see also [37]. Suppose Vn (resp. V), Cg,n (Cg) and W^'^V^C,,™ (l^^^VW^C,) 
are the Laplace exponent, the positive root (2.3) and the scale function of X^"^' {X), respectively. Because W^'^' 
is an increasing function, the measure W^'^' (dx) associated with the distribution of W^'^' (0, x] is well-defined and 
we obtain as in page 218 of [25], 

(2.23) / e-^-W^'^Xdx) = —^ . 

Because these processes are spectrally negative and tp is continuous on [0, oo), we have, by the continuity theo- 
rem, ipniP) — ^ V'(/3) for every /3 > 0. Now in view of (2.23), the convergence of the scale function holds by 
the continuity of the scale function and the continuity theorem; see [16], Theorem 2a, XIII. 1. In other words, 
Wn\x) -^ W^'''>{x) and W(^,„(x) -> W^^i^x) as n f oo for every x > 0. 

The smoothness and monotonicity properties of the scale function can be additionally used to obtain stronger 
results. The scale functions in Proposition 2.1 are in C°°(0, oo). In addition, when all the roots of Xg are different, 
its first derivative W^ is completely monotone as discussed in Corollary 2.1. 

Suppose that it is in C^(0, oo) (which holds, for example, when cr > by [12]), W^ (0+) < oo (i.e., cr > or 

n(— oo, 0) < oo) and W^' (x) < for every x > 0. Because W^ (x) > 0, we see that 

nx 



F(x):=(W^^(0+))-i / \W^^{y)\dy 



/o 
is a probability distribution and 

t-oo r o2 

/ e-^^F(dx) = (iy<^^(0+))-i -T7^ 
Jo 



+ (3W^^(0) + Wi(0+) 



/3>0. 



Therefore, noting that W^'^^ (x) = CqW^'^^ (x) + e'^'^^Wi (x) and assuming that the convergent sequence W^^ „(x) 
has the same property, we can obtain by the continuity theorem 

W^4,n(^) ^^ W/'C,(^) ^"'i Wj.^yix) ^^ W^(«)'(x), X > 0. 

The key assumption of the negativity of W^' holds, for example, for the completely monotone jump case because 
Wi is completely monotone by [33]. We can also choose the sequence Wi „ completely monotone in view of 
Corollary 2.1 because approximation can be done via hyperexponential distributions. In fact, it also means that 
W(; is C°° (0, oo) and the convergence of higher derivatives can be pursued. Even for a general jump distribution, 
the negativity of Wi' is a reasonable assumption in view of the numerical plots given by [38]. 

2.5. Obtaining overshoot and undershoot distributions. We conclude this section by giving an example how 
we can apply the explicit expression of the scale function obtained above. We consider the joint distribution of 
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overshoot and undershoot (with discounting): 



(2.24) hq{x; A, B) := E"" 



e'^^o 1 



{X _ £B,X -£A,T~<oo} 



, Be B{0,oo), AeB{-oo,0), 



focusing on the hyperexponential case addressed in Example 2. 1 . Our objective is to show an example where the 
PH (hyperexponential)-fitting can be applied; we focus on this rather simple example because we later evaluate it 
as an approximation tool by comparing the result with the simulated expectation. The results obtained here can be 
easily extended to more general fluctuation identities; see e.g. [24]. 

We assume that X is a hyperexponential Levy process and q > 0. The joint distribution (2.24) can be easily 
obtained via the resolvent measure and the compensation formula. 

Lemma 2.5. For all B e B{0, oo) and A G B{-oo, 0), 

hg{x; A,B)= f n(dn) J w'^''\x) f e-^'^'dy - f dyW^''\x -y)\, 

Jo [ JBniA+u) JBniA+u) J 

where Il{du) = XJ2T=iPj'nj^~^^^^ufor all u G (0, oo). 

After integrating with respect to 11 and assume A and B are open intervals, we obtain the following. 
Proposition 2.2. Suppose B = (b, b) and A = {—a, —a) for some < a<a and < b<b. Then 



hg{x; A, B) = XY.Pjie-'^'- " e"'^^^)/^,- ,(x; B) 



where, for each I < j < m, 



oCg^ 



Kj^q{x;B) : — 



i^'iQivj+Cg 



+ E c.. 



^i,q-^ 



e '>».'? 



iein 



1 



iVj - Ci,q) 



f,-iVj+<:q){bVx) _ ^~{r,j+Cq){b\Jx) 



{■nj-ii,g)ibAx) _ „-{rij-(,,,q)(bAx) 



g \'0 ^t-,!! 



g yij '»»,'3^ 



iVj + Cq 



By differentiating the above, we obtain the density: 



E^ 



{~X _eda,X _ eS,Tg <oo} 

^ {X . edb,X _eA,r„-<oo} 

^0 ~ ^0 



g-{ru+<:q)b _ ^'(■nj+<:g)b^ 



m 



b < X 



'(C) 'Lutein hi ' — 



V-'CC,) 



X 
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As a special case, the overshoot (resp. undershoot) density becomes, by setting B = (0, cxd) (A = (—00, 0)), 
(2.25) 

-| m 



11 



E^ 



E' 



-9^0 1 



{—X _£da,rg <oo} 
■^0 



-9^0 1 



{X _ edft.To <oo} 

^0 



■1 






6 < X, 
b > X, 



where 

K'j,g{x;{0,oo)) 



i^'iQiri.+C) 



-r,jx 



+ E ^^.9 



iein 



iVj - ihq) 



g si, 9^ g 'Ji" 



(^i + Q 



Remark 2.1. Based on our numerical experiments, it can be conjectured that 

1 1 '^ „ 1 



V''(C,)(^, + C,) 



E^v 



ieXn 



" (%■ - ^i,g) ^ 



3 



1, 



a«(i consequently 



Kj,q{x; (0, 00)) = J] Ci,g ( — 



iVj + Cq 



,m, 



^i.q"^ 



3. Scale functions for Meromorphic Levy processes 

In this section, we consider another class of spectrally negative Levy processes called meromorphic Levy pro- 
cesses. We obtain their scale functions and use these as approximation tools for a general spectrally negative Levy 
process with a completely monotone Levy measure. Similarly to the approach applied in the last section, we obtain 
the scale function using its Wiener-Hopf factorization. It has a form expressed as an countable sum of exponential 
functions which can be bounded efficiently by finite sums. 

3.1. Meromorphic Levy processes. The following is due to [22], Definition 1. 

Definition 3.1 (spectrally negative meromorphic Levy process). A spectrally negative Levy process X is said to 
be meromorphic if the following conditions hold. 

(1) The Laplace exponent '4){s) (2.1) has a countable set of real negative poles. 

(2) For every q > 0, the Cramer-Lundberg equation (2.14) has a countable set of real negative roots. 

(3) Let {rjk] k > 1} and {Ck,q] k > 1}, respectively, be the sets of the absolute values of the poles and the 
negative roots of (2.14) for fixed q > 0. Then it satisfies the following interlacing conditions: 

• • • < -Vk < -ik,q < ■■■ < -m < -C2,q < -m < -6,q < 0. 

(4) There exists a > j ^uch that rjk ~ c/c" as k ^ oo. 

(5) The Wiener-Hopf factor (2.15) is expressed as convergent infinite products 



(3.1) 



.^(^)^n^^^. 



fc=i 



Vk 
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The class of Meromorphic Levy processes complements the class of PH Levy processes described in the pre- 
vious section because it also contains those with jumps of infinite activity. As noted by Corollary 3 of [22], the 
property (3) in Definition 4.1 is equivalent to the condition that the Levy density has the form 

oo 

i=i 

This can be seen as an extension to the hyperexponential case as described in Example 2.1. For more details, see 
[22]. 

The Wiener-Hopf factor (3.1) for this process is again a rational function as in (2. 16) for the PH case. Therefore, 
this can be inverted again by partial fraction decomposition, and we have as in Corollary 1 of [22] 



(3.2) 
where 



oo 



1=1 



A 



S + Ci,c 



i,q 



"Pa (s) 



'«:<? 



6,. 



^» / zj: 1 



1 S'.g 



e»,, 



3+^ ii 



Notice by the interlacing condition that y4j ^ > for every i > 1. 



3.2. Scale functions for meromorphic Levy processes. We now obtain the scale function for the meromorphic 
Levy process. We omit the proof because it is similar to the PH case; see Appendix A.l. 

Lemma 3.1. For every q > 0, we have 

oo 

where 



i=\ 



x> 



(3.3) 






^i,q 



Q Cq + ii,q ■ 

By (2.10) and Lemma 3.1, we have, by taking the limit. 



i> 1. 



(3.4) 



Iq ■■=Y.^i,q = {^'{Cq))-' " ^C, (0) < ^- 



i=l 



The scale function can be therefore obtained by Lemma 3. 1 and (3.4). 
Proposition 3.1. For every q> {), we have 



(3.5) W^'^{x) = Y,Ci,q 



aCg^" £3 ^i,q^ 



+ W^^^(0)e^'^' = (V'(C<?))-'e^«^' - J2 Ci,ge-«''^", x > 0. 



1=1 



1=1 
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By the dominated convergence theorem thanks to the non-negativity of Ci^q and convexity of the exponential 
function, we have, for every q > and x > 0, 

oo 

1=1 
By adjusting slightly the proof of Lemma 2.3 (given in Appendix A.l), we have the following. 

Lemma 3.2. (1) The following two statements are equivalent: 

(a) cr = and n(— oo, 0) = oo, 

(b) Ei^l ^i,qii,q = ^■ 

(2) Suppose a > or n(— oo, 0) < oo. Then, for every q > 0, we have 

(3.6) - = qEA.e.,j 






where 



3.3. Approximation of the scale functions via finite sum. The scale function obtained in Proposition 3.1 is an 
countable sum of exponential functions and in reality its exact value cannot be computed. Here, we obtain bounds 
for VF(5)(-), VF(9)'(-) and Z('?)(-) in terms of finite sums. 
For every m > 1, let 

By the interlacing condition, Af^ and C^-™ are all positive and, for every i > 1, 

A^^tAi,g and cfj^ t C^g as m ^ oo. 
Now we define candidates for the upper and lower bounds of W^ respectively by 

m 
i=l 

W}^\x) := W^^\x) - 6m [e-^«^ + g-^^'+^^+i.')^] , 
for every m > 1 and x > 0, where 

m 
i=l 

which vanishes in the limit as m — )• oo by (3.4). As candidates for the upper and lower bounds of W^'^\ we also 
define 

T4^(«'")(x) := e^«^IF[^^(x) and W^^'^^'^x) := e<^''w}^\x), x>0. 

The following proposition shows that the scale functions are bounded and approximated by these functions. 
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Proposition 3.2. For every m>\ and x > 0, we have 

(3.8) W}c^\x) < W^^{x) < W^^\x) and W^^'^^^x) < W^'^Xx) < W^^'"'\x). 

Furthermore, we have 

W^^\x)^^^Wc,{x) and W^'^'^\x) ^^^W^'^\x), 
W}-^)^:,)JI^W^^{X) and W^'i^"^\x) ^^^W^'^\x), 

uniformly on x & [0, oo). 

By straightforward calculation, we can bound Z^'^) in (2.6). Let, for every m > 1, 

Z^''''")(a;) ■=l + q r W^""'''\y)dy and Z(^'™)(x) := 1 + q r W^"''''\y)dy, x > 0. 
Jo Jo 

Then by Proposition 3.2, we have Z('?'™)(x) < Z('?'™)(x) < Z^^''"^(x) and 
<Z^'^''^\x) - z!-i'"'\x) = 



W^'''^\y)-W!^i'^\y))dy 



\ _|_ c,~^m + l,qy 



dy = q5„ 



e 



m+l,(j 



.(1 _ g-im + l,qX\ 



We therefore have the following. 

Corollary 3.1 (Bounds on Z(9)). We have Z^'^'"'\x) -^ Z'^i\x) and ZJ'^^"'\x) -^ Z'-'^^x) as m ^ oo pointwise 
for every x > 0. 

We now obtain bounds for the derivative. Define, for every x > 0, 

m 



w^'^Hx) 



i=l 



w^'^Hx) := w^"'Hx) + 



max {^k,qe ^'='''^) + max {S,k,qe~ 

l<k<m k>m+l 



<k,qX) 



Om.- 



Here notice that 



max {^k,qe ^'='''''') 



j_ 

ex ' 



sm+l,g S ^ ) 



Proposition 3.3. For every m > 1, we have 

w;(™) (x) < W^^y (x) < «;(") (x) , x > 0. 
Furthermore, we have w}"^>{x) — )■ W^*^' (x) andw^"^'{x) — )■ W^'^' (x) uniformly on x > xq for any xq > 0. 

A stronger result holds when u > or n(— oo, 0) < oo. Recall in this case that < oo by Lemma 3.2-(2) and 
hence we can define 

^ m 



em:=e-^J^ei,9A-?>0, m>l 



j=i 



which vanishes in the limit as tti — ;• oo by Lemma 3.2-(2). 
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Corollary 3.2. When a > or n(— oo, 0) < oo, we have 

w;(") (x) < VF(^)' (x) < w("^) (x) A w^"^^ (x), X > 
where 

w^'^Hx) := w;(")(x) + max (efc,,e-«'=.'")5„ + e-^^+^-'^e^. 

l<k<m 
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Pj 



Vj 



1 0.029931 676.178 

2 0.093283 38.7090 

3 0.332195 4.27400 

4 0.476233 0.76100 

5 0.068340 0.24800 

6 0.000018 0.09700 



j 


Pi 


Vj 


J 


Pj 


^j 


1 


8.37E-11 


8.3E-09 


8 


0.000147 


0.0020 


2 


7.18E-10 


6.8E-08 


9 


0.001122 


0.0100 


3 


5.56E-09 


3.9E-07 


10 


0.008462 


0.0570 


4 


4.27E-08 


2.2E-06 


11 


0.059768 


0.3060 


5 


3.27E-07 


1.2E-05 


12 


0.307218 


1.5460 


6 


2.50E-06 


6.5E-05 


13 


0.533823 


6.5160 


7 


1.92E-05 


3.5E-04 


14 


0.089437 


23.304 



(i) Weibull(0.6,0.665) (ii) Pareto(1.2,5) 

Table 1 . Parameters of the hyperexponential distributions fitted to (i) Weibull(0.6, 0.665) and 
(ii) Pareto(1.2, 5) (taken from Tables 3 and 9, respectively, of [15]). 



4. NUMERICAL Examples 

In this section, we illustrate numerically the effectiveness of the PH-fitting approximation for a general spectrally 
negative Levy process. First, we use the classical hyperexponential-fitting algorithm for a completely monotone 
density function by [15] to approximate the scale function for the case with a (Brownian motion plus) compound 
Poisson process with Weibull/Pareto-distributed jumps. We approximate (2.5) and (2.24) and compare them with 
simulated results. Second, we consider, as an example of the meromorphic Levy process, the /3-family introduced 
by [20] and extend the results to the spectrally negative version of the CGMY process. 

4.1. Numerical results on hyperexponential-fitting. As noted earlier, any spectrally negative Levy process with 
a completely monotone Levy density can be approximated arbitrarily closely by hyperexponential-fitting. Here, 
we use the fitted data computed by [15] to approximate the scale function for a compound Poisson process with 
i.i.d. Weibull/Pareto-distributed jumps (with or without a Brownian motion component). Recall that the WeibuU 
distribution with parameters c and a (denoted Weibull(c, a)) is give by 

F(t) = l-e-(*/'')', t>0 

and the Pareto distribution with positive parameters a and b (denoted Pareto(a, b)) is given by 

F{t) = l-{l + bty, t>0. 

These have long-tails, namely, e''*(l — F{t)) — )• oo as t f c« for any 6 > 0. See [19] for more details about these 
distributions. Feldman and Whitt [15] constructed a recursive algorithm to approximate completely monotone 



16 



M. EGAMI AND K. YAMAZAKI 



1 







scale function 



simulation 



scale function 



simulation 



.30312 .30256 (.30013, .30498) .30977 

.46728 .46793 (.46517, .47068) .47102 

.63707 .63642 (.63370, .63914) .63862 

.81409 .81231 (.81023, .81439) .81433 



.31053 (.30797, .31309) 
.47113 (.46792, .47435) 
.63662 (.63372, .63952) 
.81546 (.81310, .81782) 



Exp(l) 





a = 1 


a = 


X 


scale function 


simulation 


scale function 


simulation 


1 


.40683 


.40420 (.401 16, .40724) 


.41734 


.41517 (.41225, .41809) 


2 


.55851 


.55473 (.55176, .55769) 


.56526 


.56394 (.56080, .56707) 


3 


.70570 


.70013 (.69767, .70259) 


.70968 


.70491 (.70229, .70753) 


4 


.85210 


.84925 (.84740, .85111) 


.85385 


.85240 (.85017, .85462) 







Weibull(0,6,0.665) 






a = 1 


a = 


X 


scale function 


simulation 


scale function 


simulation 


1 


.69953 


.70291 (.70042, .70540) 


.72042 


.72287 (.72041, .72533) 


2 


.80667 


.80782 (.80569, .80994) 


.81614 


.81729 081513, .81945) 


3 


.88181 


.88176 (.88003, .88348) 


.88631 


.88663 (.88491, .88835) 


4 


.94412 


.94364 (.94248, .94480) 


.94590 


.94747 (.94635, .94859) 



Pareto( 1.2,5) 
Table 2. Computation of E^[e''?^^^l|^->^+,,+<^|] 
The simulation column shows the mean and the 95% confidence interval. 



via scale function and simulation. 



densities in terms of hyperexponential densities. We use their results and approximate the scale functions of 
spectrally negative Levy processes with Weibull/Pareto-distributed jumps. 

We consider the Levy processes X^*^"'"") and x(p^'''='°) in the form (2.12) where Z is (i) Weibull(0.6,0.665) 
and (ii)Pareto(1.2,5), respectively. Table 1 shows the parameters of the hyperexponential distributions obtained 
by [15] fitted to (i) with m = 6 and to (ii) with m = 14. We use these parameters to construct hyperexponential 
Levy processes x(*'=i''"") and X(P^''=to) (see Example 2.1) that will be used to approximate X^*"''"") and x(p^'''='°), 
respectively. 

First, we study how accurately the scale function of x^^'^'*'"") (resp. X^^^^^^°^) can be approximated by those of 
^(weibuU) (-^(pareto)) joward tWs end, we use the identity as in (2.5) 



(4.1) 



E^ 



e-'i^bi 



{^0 >^^' 



'<oo} 



iy(^)(x)/Ty(«)(6). 



We compute the right-hand side for x^*'^'''"") (xiP"^^^")) exphcidy via (2.20) and approximate the left-hand side 
for x('^'='''"") (x(P^'''='°)) via Monte Carlo simulation based on 100, 000 sample paths. For the simulated results, 
Brownian motions are approximated by random walks with time step At = T/lOO for each interarrival time T 



ON SCALE FUNCTIONS OF SPECTRALLY NEGATIVE LEVY PROCESSES WITH PHASE-TYPE JUMPS 



17 



between jumps. In order to confirm the accuracy of the simulated results, we also calculate both sides for the 
process X^'^^p) with i.i.d. exponential jumps with parameter 1. For the scale functions of X^^^^', xi^^^'"^^^) and 
^(pareto)^ the roots ^.,g's and Q are obtained via the bisection method with error bound l.OE — 10. 

We consider the case a = and a = 1 and various values of starting point x with common parameters 
q = 0.05, A = 5, 6 = 5 and fi = 5. Table 2 gives the results. From these, we see that (4.1) for x^*'^'"'"") and 
■^{paieto) ^g approximated very precisely. Based on the accuracy for each x = 1, . . . , 4, we can also infer that 
their scale functions are approximated efficiently. Their scale functions and first derivatives for X^'^^p^ j^(weibuii) 
and X^P^"^*^'") are plotted in Figure 1. It can be easily confirmed that the results are consistent with Lemma 2.2; in 
particular, W'-'^^O) = l//i = 0.2 and P^('?)'(0+) = {q + A)/(/x2) = 0.202 for a = whereas iy(«)(0) = and 
|y('?)'(0+) = 2/(^2) = 2 for (T = 1. 

Second, we evaluate the approximation of (2.24) for x(^'='''"") (X^p^'^'")) using those for x(^'="'"") (^(pareto)^ 
Here we compute the overshoot/undershoot density 



E^ 



e-9^0 1 



{—X -£da,Tn <oo| 
^0 



and E^ 



e-i^o 1 



{X _ ed&, T„ <oo} 
^0 ~ 



These values are computed explicitly by (2.25) for x(*'='''"") (X^p^''^'")). For x^""^^''"") (x(p^'''='°)), we simulate 



E^ 



-9^0 1 



{~X _e(a-Aa/2,a+Aa/2),r„ <oo} 



/Aa and E^ 



-9^0 1 



{X _ e{b-Ab/2,b+Ab/2),T~<oo} 



/Ab 



with Aa = Ab = 0.1 by Monte Carlo simulation with 500,000 samples. These values are also computed for 

ji^(exp)_ 

Figures 2-3 show the results for the cases a = 1 and o" = with common parameters x = 5, // = 1, A = 10 and 
q = 0.05. In Figure 3, the density has a jump at the initial position x = 5 for the case a = (while it is continuous 
for the case a = I) due to the fact that x = 5 is irregular for (0, 5). As can be seen from these figures, the 
approximation accurately captures the overshoot/undershoot density for x^^'^'*'"") and X^p^^'"). The spike at the 
initial position in Figure 3 is precisely realized thanks to the closed-form expression (2.25); this would be difficult 
to realize if the scale function is approximated via numerical Laplace inversion. 

4.2. Numerical results on the /3-class and CGMY process. We now consider, as an example of meromorphic 
Levy processes, the /3-class introduced by [20]. The following definition is due to [20], Definition 4. 

Definition 4.1. A spectrally negative Levy process is said to be in the f3-class if its Levy density is in the form 



(4.2) ^(^) = C (i_g/3x)A ^{x'<o}, xeR, 

for some a>0, /3>0, c>0 and A G (0, 3). It is equivalent to saying that its Laplace exponent is 

V^(z) =[iz + ]^cj^z^ + ^ |i?(a + ^, 1 - A) - B{a, 1 - A) 

where B is the beta function, i.e., B{x, y) := r(x)r(y)/r(x + y). 

The special case cr = and /3 = 1 reduces to the class of Lamperti-stable processes, which are obtained by 
the Lamperti transformation [31] from the stable processes conditioned to stay positive; see [8, 9] and references 
therein. For the scale function of a related process, see [28]. 
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Scale function for the Exp(l) case 



Derivative for the Exp(l) case 
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Scale function for the Weibull(0.6,0.665) case 



Derivative for the Weibull(0.6,0.665) case 
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Scale function for the Pareto(1.2,5) case 



Derivative for the Pareto( 1.2,5) case 



Figure 1 . Scale functions and their derivatives. 
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Weibull(0.6,0.665) with a = l 



overshoot 

Weibull(0.6,0.665) with a = 
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Pareto( 1.2,5) with o- = 1 



" via scale function 
via simulation 
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Pareto( 1.2,5) with o- = 



Figure 2. Computation ofthe overshoot density E''[e ^^o Ir ^ _edar"<oo)]- 
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Pareto( 1.2,5) with o- = 1 
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Pareto( 1.2,5) with o- = 



Figure 3 . Computation of the undershoot density E^ [e ^^o 1 
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It can be also seen that this is a "discrete-version" of the (spectrally negative) CGMY process, whose Levy 
density is given by 

(4.3) 7r(x) = c-nyl|^.<0}, xGM. 
Indeed, if we set c = c/3^ and a = af3~^ in (4.2), we have 

(4.4) ^ (i_e/3x)A l{^<0} ^ ^^l{-<0}' ^ e M. 

See [4] for approximation of (double-sided) CGMY processes using hyperexponential distributions. 

We use the results in Section 3 to obtain the bounds on the scale functions. We consider the case q = 0.03, 
o" = 0.2, (1 = 0.1, A = 1.5, a = 3, /3 = 1 and c = 0.1 in (4.2). It is easy to verify that these have jumps of infinite 
activity (and bounded variation). First, we plot, in Figure 4, the upper and lower bounds on the scale function and 
its derivative as obtain in Proposition 3.2 and Corollary 3.2 for ?tt, = 10 and m = 100. The difference between the 
upper and lower bounds indeed converges to zero, although the convergence of the derivative in the neighborhood 
of zero is relatively slow. Second, motivated by (4.4), we take /3 to zero and see how the approximation for the 
CGMY process works. Here we set d = 3, c = 0.1 and m = 100 and use the same values as the above for the 
other parameters. Figure 5 shows the upper and lower bounds of scale function and its derivative for various values 
of /3; we can indeed observe the convergence as /3 — )• 0. 

5. Concluding Remarks 

We have studied the scale function for the spectrally negative PH Levy process and the PH-fitting approach for 
the approximation of the scale function for a general spectrally negative Levy process. Because the approximated 
scale function is given as a function in a closed form, one can analytically differentiate/integrate to obtain other 
fluctuation identities explicitly. Our numerical results suggest that the PH-fitting is a powerful alternative to the 
numerical Laplace inversion approach, particularly when the Levy density is completely monotone. 

One major challenge for the PH-fitting approach is that, without the completely monotone assumption, there 
does not exist a fitting algorithm that is guaranteed to work for arbitrary measure. The performance is certainly 
dependent on the shape of the Levy measure and it may not be a suitable approach for certain cases. However, 
there exist a variety of fitting algorithms typically developed in queueing analysis. Well-known examples are the 
moment-matching approach (e.g. MEFIT and MEDA) and the maximum-likelihood approach (e.g. MLAPH and 
EMPHT), and a thorough study of pros and cons of each fitting techniques has been conducted in, for example, 
[17, 32]. Our next step is, therefore, to apply these existing algorithms for the approximation of the scale function, 
and analyze its performance for a variety of Levy measures. Another direction for future research is the PH-fitting 
construction of scale functions from empirical data as in, for example, [4]. The closed-form expression of the 
approximated scale function can be used flexibly to identify the fluctuation of the process implied by the empirical 
data. 

Appendix A. Proofs 



22 



M. EGAMI AND K. YAMAZAKI 



" upper bound / 
" lower bound 

















upper bound -' 


- 










~ ~ ~ " lower bound 




/■ 










/ 


- 








/ 

/ 
/ 


- 














^^-^^ 


^''^ 















Bounds on VF^) (x) when m = 10 



Bounds on W^'^'^ (x) when m = 100 





















~ ~ ~ upper bound 
















~ ~ ~ " lower bound 




/ 
















/ 


250 

200 














/ 
/ 
/ 
/ 

/ 


150 
100 


s 










^/' 


/ 

/ 
/ 


50 


-- 


- =-- 


- -' " 

















- - - upper bound 
~ ~ ~ " lower bound ^ 


/ 

/ 
/ 

/ 
/ 

/ 
/ 

y 



Bounds on I^('^) (x) when ?ri = 10 



Bounds on W^'^^ (x) when ?n = 100 



Figure 4. Approximation of the scale function and its derivative for the /3-class. 



A. 1. Proof of Lemmas 2.3-2.4. By (2. 17), it is easy to verify that 
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Figure 5. Convergence to the CGMY model. 



and hence, because W^'^^ G C'^lO, oo), 
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On the other hand, different representations of (A.1)-(A.2) can be pursued. By Theorem 8.1 of [25] and (2.7), 
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for every < a < x. Its derivative with respect to a becomes 
d 
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In particular, when o = 0, the derivative with respect to x and its limit as x — )• are 
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By matching (A.2) and (A.4), Lemma 2.3 is immediate. 

For the proof of Lemma 2.4, by matching (A. 1) and (A.3), we have 
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Integrating the above and changing variables, we have 
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Lemma 2.4 is now immediate because the integral part is a lower incomplete gamma function. 

A.2. Proof of Lemma 2.5. Let N{-, ■) be the Poisson random measure for —X and X^ := mino<s<i Xg for all 
i > 0. By the compensation formula, we have 
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By using the q-resolvent kernel (see, e.g., [25]) 
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Substituting this, we obtain the result. 



dt [e-«*P^{Xi_ G S n (^ + u), Tq- > i}] = / dy 
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A.3. Proof of Proposition 2.2. Define for all B G B{<d, oo), A G S(-oo, 0) and /sT G R 

p{K;A,B):= Tl{du) e^^dy. 
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we obtain via Lemma 2.5 
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Lemma A.l. Suppose B = (6, b) and A = {—a, —a). For any K G M, we have 
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For its integral with respect to IT, 

n(dn)e^«-^+")^^) = r n(du)e^(-^+") + f n(du)e^^ = e^^ rTl{du + b)+ f n(du + a)e^", 

J a+b J a+b Ja Jb 



J a+b J a+b J a+b 

ra+b ra+b /"a+b 



I a+b 



n(dn)e^«-'^+^)^^) = r "n(du)e^^+ T n(du)e^(-"+") = e^** P U{du + b) + f U{du + a 

J a+b J a+b Ja Jb 



^^^Ku_ 



26 M. EGAMI AND K. YAMAZAKI 

Therefore 



fb ca fb 
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We are now ready to prove Proposidon 2.2. By Lemma A.l, 
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^.4. Proof of Proposition 3.2. Notice, for every m > 1, that 
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Therefore we have the bounds for W(^^ in (3.8). The bounds for W^'^^ are immediate by muhiplying e''"'^. Finally, 
the convergence results hold because dm ^ ^ and e~'»''^ + e~(^9+^™+i''?)^ and 1 + e~'^™+i'J^ are bounded uniformly 
in X > 0. 

A.5. Proof of Proposition 3.3. The] 
we have by (A.5) 
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which shows the first claim. For the 
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A.6. Proof of Corollary 3.2. For every x > 0, 
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where the last inequaUty holds because 

^1 Sr^ c. A. -a ^iSr-c. A. ^a ^iSr-c. /iM 



=m+l ^ k=l 

This together with Proposition 3.3 shows the claim. 
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